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We generalize our recent approach to reconstruction of phase dynamics of coupled oscillators from data [B. 
Kralemann et al., Phys. Rev. E, 77, 066205 (2008)] to cover the case of small networks of coupled periodic 
units. Starting from the multivariate time series, we first reconstruct genuine phases and then obtain the 
coupling functions in terms of these phases. The partial norms of these coupling functions quantify directed 
coupling between oscillators. We illustrate the method by different network motifs for three coupled oscillators 
and for random networks of five and nine units. We also discuss nonlinear effects in coupling. 

PACS numbers: 05.45.Xt Synchronization; coupled oscillators 
05.45.Tp Time series analysis 



Many natural and technological systems can be 
described as networks of coupled oscillators. A 
typical problem in their analysis is to find dynam- 
ical features, e.g., synchronization, transition to 
chaos, etc., in dependence on the properties of 
oscillators and of the coupling. Here we address 
the inverse problem: how to find the properties of 
the coupling from the observed dynamics of the 
oscillators. This may be relevant for many exper- 
imental situations, where the equations of under- 
lying dynamics are not known. We present here 
a method which is based on invariant reconstruc- 
tion of phase dynamics equations from multivari- 
ate observations (at least one scalar oscillating 
observable of each oscillator must be available) 
and includes several algorithmic steps which are 
rather easy to implement numerically. We illus- 
trate the method by numerical examples of small 
networks of van der Pol oscillators. 



I. INTRODUCTION 

Dynamics of coupled self-sustained oscillators attracts 
vast interest of researchers, both due to a variety of non- 
trivial effects and to numerous applications in physics, 
chemistry, and life sciences. In modeling, one usually 
focuses on the relations of the dynamical features, such 
as synchronization patterns, to the coupling structure of 
the underlying network. In experimental data analysis, 
an inverse problem typically arises, where one wants to 
reveal the interactions from the observed dynamics. A 
simple consideration demonstrates that this task is re- 
ally challenging and highly nontrivial: indeed, computa- 
tion of correlations and or synchronization index between 
the oscillators does not solve the problem, because the 
interactions are generally non-reciprocal while both the 
correlation and the synchronization index are symmetric 
measures. 



There are exist two approaches for the determination 
of directional coupling. The first one exploits various 
information-theoretical techniques and quantifies an in- 
teraction with the help of such asymmetric measures 
as transfer entropy, Granger causality, predictability, 
eteii~— . The second approach, developed in our previ- 
ous publications^ 7 -, see also£, is based on estimation of 
phases from oscillatory time series and subsequent re- 
construction of phase dynamics equations and quantifica- 
tion of coupling by inspecting and quantifying the struc- 
ture of these equations. The disadvantage of the second, 
dynamical, approach in comparison to the information- 
theoretical one is that it is restricted to the case of weakly 
coupled oscillatory systems. However, for this widely en- 
countered case the dynamical approach has clear advan- 
tages: it works with phases, which are most sensitive to 
interaction, and yields description of connectivity which 
admits a simple interpretation. This consideration is sup- 
ported by several comparative studies^ '. (Notice that 
there is an intermediate group of techniques which ap- 
ply informatical measures to time series of phase a 7 ' 11 " — .) 
The dynamical approach has been tested in a physical 
experiment^ and used for analysis of connectivity in a 
network of two or several oscillators in the context of 
cardio-respiratory interaction 15 , brain activity.^—, and 
climate dynamics^. 

Recently, we have essentially improved the dynamical 
approach by suggesting a technique for invariant recon- 
struction of the phase dynamics equations from bivariate 
time series^—. Invariance in this context means that 
reconstructed equations do not depend on the observ- 
ables used (at least for a wide class of observables) . This 
technique was tested on numerical examples as well as in 
physical experiments with coupled metronome o 20 ' 21 and 
with electro-chemical oscillators 2 ^. In this paper we ex- 
tend the technique to cover the case of small networks of 
interacting oscillators. The paper is organized as follows. 
In Section |TT] we discuss phase description of oscillator 
network. In Section ITTT1 we describe the basic method; it 
is illustrated in Sections IIVI and W\ We summarize and 
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discuss our results in Section fVTl 

II. PHASE DESCRIPTION OF OSCILLATOR NETWORK 

We assume that an individual oscillator is described 
by variables x which satisfy a system of ODEs x = G(x) 
and that this equation system possesses a stable limit cy- 
cle. The latter can be parameterized by phase tp which 
grows uniformly in time, <p = u — const, where ui is the 
oscillation frequency. Notably, this phase parameteriza- 
tion is unique (up to trivial shifts) and invariant with 
respect to variable transformations^ 4 ^. 

A network of N oscillators can be represented as 

x fe = G fc (x fc ) +H fe (xi,x 2 ,...) , k = l,...,N, (1) 

where coupling functions generally depend on the 
states of all oscillators. If variables Xj are absent in 
we say that there is no direct coupling from I to k. Note 
that generally the oscillators can be quite different, e.g., 
even the dimensions of state variables X& can differ. Of 
course, the functions can be different as well. 

Assuming weakness of the coupling, one can reduce 
network equations Eq. fl} to phase equations 

ip k = U) k + h k ((pi,ip 2 ,...) , k = l,...,N , (2) 

using a perturbation technique (se o 24 ' 25 for details). 
Then, in the first approximation, h k depends explicitly 
on phase pi only if Hfe depends explicitly on x/, i.e. if 
there exists a direct link / — > k. However, in the higher 
approximations, the r.h.s. of Eq. ^ may contain high- 
order terms which depend not only on phases of directly 
coupled oscillators, but on the other phases as well. Thus, 
generally speaking, phase oscillators Eq. @ are all-to-all 
coupled. Nevertheless, as shown in Section IIIII below, 
we can decompose the coupling functions hk into several 
components, and comparing the norms of these compo- 
nents, characterize the partial couplings between partic- 
ular nodes. In fact, we attribute links only to those nodes 
for which the corresponding norms are sufficiently large. 
Reconstructed in this way, the structure of the phase net- 
work ^ correctly (although not exactly) represents the 
original network (fTJ). 

To specify the relation between the original network 
structure and that of the phase oscillators, we introduce 
the following terminology. If the function Hj, can be writ- 
ten as Hfc = Ylj^k Hkj(x-k,Xj), we call this type of cou- 
pling pairwise. Other terms, containing a combination 
of several variables acting on fc, e.g. ifjyj(xfc,Xj,Xj), we 
call cross-coupling terms. If the coupling is pairwise and 
weak, then a perturbative phase reduction yields, in the 
leading order, only pairwise terms like h k j(ipki <Pj) m t ne 
system Q . As discussed above, cross-coupling terms in 
the phase system ^ which depend on two or more driv- 
ing phases, appear either due to cross-coupling terms in 
the original system ([1} or, if the latter has pairwise cou- 
plings only, due to high-order combinational terms which 
appear in the perturbative reduction to phase dynamics. 



We emphasize that validity of the phase equations ([2]) 
is not restricted to the case of very small coupling. In- 
deed, these equations are valid as long as an attracting 
invariant torus exists in the phase space of the full sys- 
tem {!]) . This torus can also exist for parameters of cou- 
pling when a perturbation approach is not valid anymore. 
In this case, Eq. <j2j) cannot be derived from Eq. (TTJ) by 
means of a perturbation technique, but can be, however, 
reconstructed numerically from multivariate time series. 
Here, of course, the relation between networks ([2]) and 
([1]) is nontrivial; in particular, some nodes may appear 
as connected in J5]) while being disconnected in (fT|). 

The main idea of our approach is to reconstruct the 
phase dynamics ((2|) from multivariate time series; here it 
is assumed that fcth component of the series represents 
the output of the fcth oscillator. Before proceeding with 
the description of the method, we make a remark on the 
range of its validity. An important dynamical regime in 
system ^ is that of full or partial synchrony, when the 
dynamics of the phases reduces to a stable torus of lower 
dimension (partial synchrony) or to a stable limit cycle 
(if all oscillators are synchronized) . In case of synchrony 
our technique fails, because the data points do not cover 
the original ./V-dimensional torus and the coupling func- 
tions cannot be reconstructed. Sufficiently strong noise 
could cause deviation of the trajectory from the synchro- 
nization manifold and help to infer the dynamics, but 
discussion of this case goes beyond the scope of this pa- 
per. 



III. RECONSTRUCTION OF PHASE DYNAMICS FROM 
DATA 

The main assumption behind our approach is that 
we have multivariate observations of coupled oscilla- 
tors, where at least one scalar oscillating time series 
Uk{t) — 2/fc(xfc(t)) is available for each oscillator. The 
first step is to transform this observable into a cyclic ob- 
servable. Typically this is done via construction of a two- 
dimensional embedding (y k , y k ), where y can be, e.g., the 
Hilbert transform of y(t), se o 21 ' 25 for a detailed discus- 
sion. Alternatively, one can use for y the time derivative 
of y. (In all cases the data usually requires some pre- 
processing, e.g., filtering). If the trajectory in the plane 
(yk,lik) rotates around some center, one can compute a 
cyclic variable 9i(t), e.g., by means of the arctan func- 
tion. As has been in detail argued in2I, in this way we 
obtain not the genuine phase pk of the oscillator which 
enters Eq. but a cyclic variable, or protophase. The 
latter depends on the particular embedding and on the 
parameterization of rotations, and is therefore not in- 
variant. This can be seen already from the analysis of an 
autonomous oscillator: the described procedure generally 
yields a protophase 9(t) which rotates not uniformely, but 
obeys 

o = f(0) ■ 
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The transformation to the uniformly rotating genuine 
phase ip reads 



tp{0) = LO 



d9' 

W) 



where uj = 2n 



is the frequency of os- 



cillations. For practical reason it is convenient to deter- 
mine not the function /, but the probability density of 
the protophase a(9) — ^jjm ■ Thus, the transformation 
from the protophase to the phase ip can be written as 



ip(0) = 2tt / dO' a{9') . 
Jo 

Hence, determination of the genuine phase from the data 
reduces to a problem of finding the probability distri- 
bution density of the obtained protophase 9, which is a 
standard task in the statistical data analysis. Practically, 
one can use either a Fourier representation of the density 
a(9), as in^l, or a kernel function representation of a(9). 

After the transformation 9k(t) — )• <Pk{t) is performed, 
we reconstruct the phase dynamics making use of the 
fact that the time derivatives tpk are 27r-periodic func- 
tions of the phases, in accordance with Eq. @. These 
functions can be obtained from the observed multivariate 
time series of phases $^ (t) with the help of a spectral rep- 
resentation technique^ or by means of a kernel function 
estimation. Here we present the formulae for the spectral 
technique (see Sec. IV-A ir>2i for details), generalized for 
the case of more than two oscillators: 



ff>i,¥> 2 ,...) 

F u ,d = ^2 flut]- ex P(^l</>i + il 2f2 +•••)) 
h,h,... 

1 f^k(T) 

= tJ o d$fc ex p(^i $ i + ^2 +...) 
fil,.. = hf At ex p^ i$i + +•••)• 

As a result, we obtain the system of the phase dynamics 
equations in the form 

^ = fW( w ,%...)= J2 ^L..exp(zZm+^ 2 +. 

h,h,... 

(4) 

The coupling functions F^ k > describe all types of inter- 
actions: pairwise, triple-wise, quadruple- wise, and so on. 
To characterize them separately, we calculate the norms 
of different coupling terms (the partial norms), as fol- 
lows. The term J 7 ^ is a constant (phase- independent) 
one, it corresponds to the natural frequency of oscilla- 
tions. The pairwise action of oscillator j on oscillator k 
is determined by those components of F^ which depend 
on phases ipk and (pj only. We quantify this action by the 



partial norm ■N'uj ', note that here the upper index corre- 
sponds to the order of interaction (pairwise here). The 
partial norm can be computed as 



(fc) 



0,...,0,ife,0,... 1 0,Zj,0,...) 



(5) 



Correspondingly, the joint action of oscillators j, m on os- 
cillator k is determined by the cross-coupling terms con- 
taining three phases pki <p>j, <p m - This action is quantified 
by the following partial norm: 



E 



(h) 



^0 



. 1 0,/fc,0,...,0,/j,0,...,0,i m ,0,...) 



Similarly we can compute the partial norm 



(6) 
and 



so on. 

(k) 

Now we discuss the terms TJ'^ which depend 
only on one phase pk and therefore cannot be attributed 
to any interaction. In Refill we used an additional trans- 
formation ip (p' which eliminated these terms, using 
some additional assumptions on the structure of cou- 
pling. In examples presented below we checked that the 
difference between the formulations in terms of tp and <p' 
is rather small. Therefore we use the technique described 



above and neglect the small terms J-, 



O) 

0,...,0,i fc ,0,...,0 - 



IV. CASE STUDY: THREE COUPLED VAN DER POL 
OSCILLATORS 

In this section we test the presented technique on a 
system of three coupled van der Pol oscillators: 

Hi - - x\)x\ + Lofxi = A [0-12(0:2 + x 2 ) + <7\ 3 {x 3 + x 3 )} 

+ o\\r\xix z , 

x 2 - fi(l - x\)x 2 + u\x 2 = A [(721 (£1 + x\) + 0-23(2:3 + aj 3 )] 

+ V22ll x lX3 , 

x 3 - /i(l - £3)2:3 + W3X3 = A [0-31(0:1 + x\) + a 32 (x 2 + x 2 )} 

(7) 

Here the matrix <Jki, with entries zero and one, deter- 
mines the coupling structure (topology) of the network: 
for non-diagonal terms, au = 1, if there is forcing from 
oscillator I to oscillator k, and <7jy = 0, otherwise. Di- 
agonal terms Okk determine the presence or absence of 
cross-coupling. Parameters A and r\ describe the inten- 
sity of pairwise and cross-coupling, respectively. In all 
examples of this section we use the values of parameters 
fi = 0.5, ui = 1, w 2 = 1.3247, uj 3 = 1.75483. Equa- 
tions ([7]) have been solved by the Runge-Kutta method, 
and the time series of Xk , ik were used to calculate the 
protophases as 9k = — arctan(ife, Xk)- The data sets con- 
sisted of 10 6 or 10 7 points sampled with 0.01 time step. 
The analyzed coupling configurations are schematically 
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FIG. 1. Different configurations in a network of 3 oscilla- 
tors. (a)-(h): pairwise interaction. Panel (i) illustrates cross- 
coupling: an arrow from the center to, say, first oscillator 
reflects the term ~ £22:3 with an = 1 in Eq. (0, etc. 



given in Fig. [TJ Note that we chose the oscillator frequen- 
cies and parameters of coupling so that network remains 
asynchronous. 

In the first set of tests we considered only pairwise cou- 
pling in Eqs. (0, i.e. we took <7kk = 0. The results for 
A = 0.05 and A = 0.15 are presented in Figs. l2|3| respec- 
tively. Here in the left column we schematically show the 
reconstructed coupling configuration. The corresponding 
partial norms are coded by the width of arrows which link 
the nodes. Only norms which are larger than 10% of the 
maximal norm are shown here. The numerical values 
of all reconstructed norms can be seen in the tables in 
the right column. Here we show the coupling norms for 
the pairwise coupling ([5]) as non-diagonal terms, and the 
triple-phase coupling according to © as diagonal ones. 
The coupling terms, which are present in the original sys- 
tem (J7J) (i.e. those with a^i = 1) are shown in boxes. We 
see that in all cases these terms are definitely larger than 
other entries of the reconstructed coupling matrix, i.e. 
the reconstruction works pretty well in all cases. Com- 
paring the results for smaller and larger coupling (Fig. [2] 
and respectively) , we see that larger coupling strength 
leads to appearance of cross-coupling terms for phases 
that are not present in the original system (|7|). 

In the second set of tests we included the cross- 
couplings by setting Ukk = 1 and 77 = 0.1. So, each 
configuration shown in Fig. [TJa-h) was complimented by 
the links according to Fig.QJi). The results are presented 
in Fig. |U We see that basically the coupling structure 
is correctly reproduced by the method, although due to 
a joint action of different terms some couplings, not pre- 
sented in the original system (J7J, do not fall below 10% 
of the largest norm. 

Summarizing the results of this section, we conclude 
that the connectivity of a weakly and pairwisely coupled 
network can be correctly revealed. For stronger coupling, 
all really existing in the original network links are cor- 
rectly revealed and some additional links appear. How- 
ever, the latter are not necessarily spurios. Indeed, as 
discussed above, not very weak pairwise coupling in the 
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FIG. 2. Reconsructed network configurations, to be compared 
with the true configurations, shown in Fig.fjja-h). Each panel 
shows the partial norms of reconstructed coupling functions 
as a table; the numbers in boxes correspond to truly existing 
links, so ideally the numbers without the boxes should be 
zero. Pictures present the largest norms only (see text); the 
values of the norms are coded by width of linking arrows. The 
arrows from the center to the first oscillator (see panels (b), 
(f), and (h)) reflect the joint action from oscillators 2 and 3 
which cannot be decomposed into pairwise actions 2 — ► 1 and 
3 — > 1. Coupling strength A = 0.05 and r\ = 0. 



original network may lead to additional connections in 
the network of phase oscillators and to cross-coupling 
terms. Our results are in full agreement with this ar- 
gument. 



V. CASE STUDY: 
OSCILLATORS 



NETWORKS OF FIVE AND NINE 



In this section we report on the results for small net- 
works of N = 5 and N — 9 coupled van der Pol oscilla- 
tors. We simulated the systems with pairwise interaction 
only: 

Xk—n(l—xl)xk+UiXk = AyVfcifa cosafci+i; sina fci ) . 

l^k 

(8) 

As above, fi = 0.5, A is intensity of coupling, and au is 
the Nx N connectivity matrix (aki = 1 if oscillator I acts 
on oscillator k, and <7ki = otherwise). An additional 
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FIG. 3. Same as Fig. H but for A = 0.15. 
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parameter au describes the phase shift in the coupling, 
so that generally the latter can be both attracting or 
repelling. 

We note that computational efforts and data require- 
ments for reconstruction of a network of TV oscillators 
grow rapidly with N, so that full reconstruction for 
N > 3, though theoretically possible, see Eqs. (|3I4|) . be- 
comes practically unfeasible. On the other hand, in the 
previous section we have verified that weak pairwise cou- 
plings can be reliably recovered. Therefore, for networks 
described by Eqs. ([8]) we reconstruct pairwise coupling 
functions only, for all pairs of network units. 

We performed a statistical analysis of the model ([5]): 
for many randomly chosen values of frequencies u)k, pa- 
rameters a.ki , and connectivity matrices au we analysed 
pairwise interactions in the networks. For N = 5, the 
frequencies ujk have been chosen uniformly distributed 
in (0.5,1.5), and the number of incoming links for each 
node was 2 (i.e. J^i a ki — 2), while the links have been 
chosen randomly (the incoming and outgoing links were 
chosen independently). For N — 9, the frequencies uj 
were distributed uniformly in (0.5,2.5), and the number 
of incoming links was 3. The intensity of coupling in both 
cases was A = 0.15, the phase shift ctki was distributed 
uniformly in [0, 27r). 

For each network we computed N protophases accord- 
ing to 6k — — arctan(ife, WfcXfc) and then transformed 
them to phases ifk- Next, synchronization analysis have 
been performed: for all pairs we calculated the synchro- 



FIG. 4. Same as Fig. [21 but for combination of pairwise and 
triple-coupling, A = 0.1, r\ — 0.1. Panels (a-h) correspond to 
combination of pairwise coupling as shown in Fig. [Ha-h) with 
the triple-coupling as in Pig- [Hi) ■ Panel (i) shows the results 
for purely triple-coupling. 



nization index \{e i( - tpk ~ lfil ' ) )\. Since our technique does not 
work in case of synchrony, we excluded from the further 
consideration all networks where synchronization index 
was larger than 0.5 at least for one link. 

For each obtained non-synchronous network we fitted 
the time derivatives (fi k by J2k^i F ki(<Pk,<Pi), where F kl 
depend on two phases only. Next, we computed norms 
Afki of all Fki; these norms represent reconstructed con- 
nectivity of the network. Performing the analysis with 
a large ensemble of non-synchronous networks, we sepa- 
rated Afki into two classes: first class contains those val- 
ues of Afki for which cry = 1, i.e. the connections between 
nodes k, I in the network ([5]) really exists. Otherwise, if 
<Jki = 0, the value A4z belongs to the second class. The 
obtained results are summarized in Fig. [5l Ideally, one 
would expect that all norms in the first class are much 
larger than those in the second one (cf. Figs. l2l3[ ). We see 
that there is a clear, although not ideal, separation be- 
tween connected and non-connected links, what allows us 
to conclude that generally our technique correctly repro- 
duces the network structure. However, a more detailed 
statistical analysis is needed before a "blind" application 
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analysis for a given network remains an open problem. 
A possible solution might be reconstruction of phase dy- 
namics for several or all triplets of oscillators and com- 
parison of the terms, dependent on three phases, with 
those dependent only on two. If the triple terms are 
much smaller in norm than pairwise ones, than the pair- 
wise analysis may be sufficient. 



FIG. 5. Probability distribution densities p(A/") of the norms 
of reconstructed couplings: black solid line: really present 
connections [o^i = 1), dashed red line: absent connections 
(<Tfcj = 0), for networks Q of 5 (a) and 9 (b) oscillators. 



of the method would be possible. 



VI. DISCUSSION AND CONCLUSION 

In summary, we have presented a technique for recon- 
struction of the phase model of a network of coupled limit 
cycle oscillators. From the reconstructed model we in- 
fer directional couplings of the network. We have tested 
the technique on small networks of van der Pol oscilla- 
tors with pairwise and cross-coupling. We have shown, 
that presence and direction of existing links can be reli- 
ably revealed. However, if the coupling is not very weak, 
the technique yields additional links which are absent 
in the original network. The differentiation between the 
truly existing and additional links is not perfect. Now 
we discuss where the additional links come from. First, 
some additional links may be numerical artefacts (spuri- 
ous links). They may appear due to noise, closeness to 
sychrony, neglection of the amplitude information, etc. 
Second, some links, absent in the model in terms of state 
variables, should indeed be present in the phase model 
and are, therefore, not spurious. If we asume that a 
matrix au determines the structural connectivity of the 
original network model ([1]), then the phase model ([2]) for 
this network would have an effective connectivity matrix 
Sfe; with a larger number of links and possibly with many 
cross-couplings. It is this effective connectivity that is 
determined in our approach, not the original one. More- 
over, here one should distinguish between the structural 
connectivity (described by matrices a, E) and the func- 
tional one, which is described by similarities in the dy- 
namics (e.g. via correlations or synchronization). The 
functional connectivity (this term, typically used in the 
context of analysis of brain activity, is widely understood 
as a correlated time behavior) results from the dynam- 
ics, and may only loosely correspond to the structural 
one» 

Finally, we discuss the perspectives of the analysis of 
networks of many oscillators. As we have demonstrated, 
generally a pairwise analysis yields good results only if 
the connections are pairwise and weak. In the case of 
cross-coupling (see e.g., Fig. [TJi) ) , pairwise analysis does 
not help. Determination of applicability of a pairwise 
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